Modeling pulsativity in the hypothalamic–pituitary–adrenal hormonal axis

A new mathematical model for biological rhythms in the hypothalamic–pituitary–adrenal (HPA) axis is proposed. This model takes the form of a system of impulsive time-delay differential equations which include pulsatile release of adrenocorticotropin (ACTH) by the pituitary gland and a time delay for the release of glucocorticoid hormones by the adrenal gland. Numerical simulations demonstrate that the model’s response to periodic and circadian inputs from the hypothalamus are consistent with those generated by recent models which do not include a pulsatile pituitary. In contrast the oscillatory phenomena generated by the impulsive delay equation mode occur even if the time delay is zero. The observation that the time delay merely introduces a small phase shift suggesting that the effects of the adrenal gland are “downstream” to the origin of pulsativity. In addition, the model accounts for the occurrence of ultradian oscillations in an isolated pituitary gland. These observations suggest that principles of pulse modulated control, familiar to control engineers, may have an increasing role to play in understanding the HPA axis.


Scientific Reports
| (2022) 12:8480 | https://doi.org/10.1038/s41598-022-12513-w www.nature.com/scientificreports/ Both of the hypothalamo-pituitary hormones CRH and ACTH are released as pulses 30,31 . Corticoid hormones are released more smoothly, but inherit a pulsatile pattern from CRH and ACTH 31 . For a number of years it was assumed that the center of ultradian pulsativity of HPA was the hypothalamic nucleus (see e.g. 31,32 ). Although this hypothesis was confirmed for the reproductive and growth hormones, it was disproved for the HPA axis (see a discussion in 33 ). It was shown in vivo that an ultradian CORT rhythm exists, even when the hypothalamus is surgically disconnected from the pituitary 34 . Another reason to doubt the leading role of the hypothalamic frequency is that the frequency of the CRH release was found to be three times higher than that of ACTH 30,35 . Finally it has been shown that ultradian oscillations in CORT persist in the presence of a constant level of CRH 36 . Taken together these observations strongly indicate that the pulse generator in the HPA axis is sub-hypothalamic.
Our discussion is organized as follows. First, we build our model upon realistic mathematical models developed previously for the HPA axis which emphasize that pulsatile glucocorticoid production arises due to a subhypothalamic pulse generator as a result of the interplay between the pituitary and the adrenal gland 8,11,16,20,[36][37][38][39] . An integrate-and-fire mechanism is used to illustrate the effects of pulsativity of the anterior pituitary gland. Then, using numerical simulations, we demonstrate that this model reproduces the experimentally observed patterns of ultradian and circadian oscillations as the hypothalamic input to the pituitary gland is varied. Next we show that, in contrast to previous models 8, 16 , the time delay is not critical for oscillatory behavior, but merely introduces a phase shift. Finally, we show how ultradian oscillations can arise in an isolated pituitary gland, i.e. when the pituitary is disconnected from both of the hypothalamus and the adrenal gland 34 . Notice that if the level of CORT is zero, the pituitary-adrenal pulse generator proposed in 16 does not induce periodic oscillations of ACTH, which were observed in vitro in an isolated human pituitary 34 . Despite the simplicity of our model, it not only captures the salient features of the ultradian and circadian rhythms generated by the HPA axis, but also explains previously unexplained observations.

Model
Background. Our model builds on a realistic mathematical model for the HPA ultradian pulsativity proposed in 8,11,16,20,[36][37][38][39] . Following Walker et al. 16 we assume: (1) the CRH level is constant during the system's transient time; (2) the impact of glucocorticoids on the CRH production is not important for the system's dynamics and can be neglected; (3) a discrete time delay was introduced in the term describing the CORT production. It was shown that for some physiologically reasonable values of the CRH constant level and of the delay this modified 3D system exhibited periodic oscillations that could be interpreted as ultradian; (4) When the hypothalamic input is periodic, a realistic circadian rhythm is obtained.
The mathematical model that we propose has the following significant distinctions from the previously known models: (1) describing the corticotroph pulse generator with the help of an electrical action potentials was previously proposed in 40,41 . However these works used a Hodgkin-Huxley-like formalism, which is much more complicated and employs more equations than the integrate-and-fire formalism that we use; (2) the models based on 16,22 contain a nonlinear equation that describes dynamics of glucocorticoid receptors (GR) in the pituitary. Our model is more parsimonious and does not consider GR; (3) in case the inputs from the adrenal gland and from the hypothalamus to the pituitary are lacking, the concentration of ACTH vanishes in the GR-oriented models. In our model oscillations of ACTH are supported even for zero inputs; (4) in GR-oriented models oscillativity of the hormone concentrations is attained by the introduction of a time delay into the equation describing the CORT release. Our model remains oscillative without such delay.
Additional assumptions. We make the following additional assumptions: (1) The center of HPA ultradian pulsativity is located in pituitary 8,11,16,36,42 ; (2) ACTH and CORT are released in pulses with an ultradian frequency of approximately one pulse per hour 8 ; (3) The pituitary hormonal cells generate some membrane electrical potential, just like neurons, and this potential controls the amount of released ACTH 43,44 ; (4) The pre-synthesized ACTH is accumulated in secretory vesicles near the cell membrane. It is released very rapidly in response to hypothalamic stimulation 8 ; (5) When ACTH reaches the adrenal gland, it launches the process of CORT synthesis. Thus, unlike ACTH, CORT cannot be released immediately, but requires some time for synthesis 9 . Following 8 , we took this circumstance into account by choosing the time delay of 15 min in the feedforward connection between the pituitary and the adrenal gland; (6) The level of CORT follows the circadian profile. It is high during the active time and low in the sleep period. A disruption of the circadian rhythm nevertheless preserves the ultradian rhythm 34,45 ; (7) The release of CRH is also pulsatile with an approximate frequency of three pulses per hour. While the CRH pulse frequency is rather steady during the diurnal period, the amplitudes of CRH pulses vary significantly, following the circadian rhythm 30,42 ; (8) The increase of the constant CRH level implies an increase of the ultradian frequency 8 . The blockade of endogenous CRH results in a significant reduction of the amplitudes of ultradian pulses 42 . (9) CORT negative feedback at the hypothalamus is not Figure 1. Schematic representation of the HPA axis studied in this paper. Arrows and bar-headed lines indicate excitatory and inhibitory connections, respectively. Following 16  Mathematical formulation. Following the assumptions of the previous section, we suggest the following mathematical model. Consider a system of differential equations with delay. The two main variables involved are A(t) and C(t) that represent concentration levels of ACTH and CORT in blood, respectively.
(1) Equation (1) relates to the adrenal gland and follows 16 . The concentration C(t) satisfies a differential equation with delay where k c is a degradation rate of CORT and the coefficient k ca relates to the driving (feedforward) signal from the pituitary to the adrenal gland. The delay τ reflects the time required for the CORT synthesis. (2) Equations (2)-(5) describe the dynamics of the ACTH release by the corticotropic cells in the anterior pituitary. The spiking dynamics of these cells are quite complex and include spiking, bursting and irregular spiking patterns. These complexities arise because of the effects of BK (big potassium) channels on the dynamics of conductance-based models (see Discussion). Here we make the simplifying assumption that the pituitary release of ACTH can be described by an integrate and fire mechanism. The pulse generator located in the anterior pituitary is modeled with two equations, the first of them is ordinary differential and the second one is a differential equation with jumps (impulses): The variable V(t) is a charging-discharging membrane potential related to pituitary cells (cf. 44 ). The pulsation times are defined from the recursion where is a given firing threshold. After the firing the potential V(t) resets to zero: Here V (t + n ) is the right-sided limit of V(t) at the point t n . Together (3)-(5) implement an integrate-and-fire scheme with a leakage coefficient k v . a threshold and an excitatory input (cf. 44 ). The greater is I(t), the faster V(t) reaches the threshold, so the firing frequency increases. If the threshold is never reached, there will be no impulses (the system approaches an equilibrium). All the jumps of V(t) are of the same value (equal to ).
The function F a (C) describes a feedback from CORT to ACTH and V. It is decreasing in C to ensure an inhibition of the ultradian amplitude and frequency when the level of CORT increases. With the increase of C(t) the right-hand side of (2) and the excitatory input I(t) decrease, so the level of A(t) also decreases (is inhibited) and the impulses of V(t) become sparser. In our simulations the function F a (C) is taken to represent repression Michaelis-Menten kinetics where F 0 , h c are some positive parameters.
Notice that the presence of feedback F a (C) is not critical for oscillativity of ACTH. Namely, if we neglect the inhibitory effect of CORT by setting F a (C) ≡ F 0 = const , the function A(t) still may oscillate. The main role of this feedback is an additional (to the hypothalamic input) modulation of amplitudes and frequencies.
The non-linearity F v (V ) is empirical, it relates to the shape of a single ACTH pulse. These pulses are rapidly released in a train, following exactly the frequency determined by V(t). In our simulation experiments we took the following function used in neuroscience for an action potential: either the exponentially decreasing function F v (V ) = e −k s V (see 10,47 ), or the alpha function F v (V ) = V e −k s V (see 48 ), where k s is a positive parameter determining the decay rate.
The parameter k a is a degradation rate of ACTH and can be obtained from experimental data. The rest of parameters k ah , k vh , A 0 , V 0 are all positive, they are adjusted empirically to fit the known HPA dynamics.
The structure of Eqs. (1)-(5) seems to be new and develops the mathematical scheme put forward in 5 . Unlike the mathematical formalism of 16 , our Eqs. (1)-(5) are not just differential, but functional-differential and contain impulses (cf. 49 ).
(3) The function H(t) is a synergetic input from the hypothalamic nuclei that combines the CRH drive H r (t) (from PVN) and the exogenous circadian signal H c (t) (from SCN).
Following 16 , we neglect a feedback from CORT to the hypothalamus. Thus H r (t) is defined periodic with a given period T r and a given degradation rate k r : Here ω = 2π/1440 is the circadian frequency and t 0 shifts the peak of H c (t) to the time of the expected maximum of circadian activity (cf. 19 ).
The period T r of H r (t) is usually about an hour, so H r (t) can be considered high-frequency when compared with the slow circadian drive H c (t) . At a short interval of simulation (of several hours) the circadian input can be considered constant, and (if H(t) is properly normalized) we can take H(t) = H r (t) neglecting circadian variations. In a more general situation we consider H(t) as a product By this multiplication the two signals are superimposed and the high-frequency signal H r (t) modulates the low-frequency circadian drive to obtain a cumulative drive H(t) of a more complex spectrum.
In system (1)-(9) the coordinates A(t), C(t) are continuous in time, and V(t) has jumps at the points t n defined from functional relationships (4), (5). Some mathematical properties of the system are described in the Supplementary Material. Parameters used for simulation. For definitness, assume that A(t) and C(t) are measured in pg/ml and ng/ml, respectively, and time is measured in minutes. The parameter values used for computer simulations are presented in Table 1.
The first five parameters in Table 1 are taken from 36,50 . Namely, from the supplementary material to 36 , the ranges for the half-life times ( T 1/2 ) of ACTH and CORT can be taken 0.5-1 min and 7.2-10 min, respectively. (The authors refer to experiments in 23,51 .) The decay rate of an exponentially decreasing signal can be calculated as ln(2)/T 1/2 . This gives the ranges 0.693-1.386 min −1 for k a and 0.069-0.96 min −1 for k c . Following the main text of 36 , we take the delay of τ = 15 min. In 50 the half-life time of CRH is defined as 25.3 ± 1 min, that gives the decay rate in the range 0.0264-0.0285 min −1 . (Notice that experimentally found values of hormonal halflives vary significantly from one publication to another, so this is one of the possible choices.) The rest of the parameters (such as coupling coefficients and gains) were adjusted to achieve the best dynamic fit to the known activity of the HPA axis (cf. 19 ).
Notice that the parameters that fit a real date may significantly vary from one human individual to another. This is even more true for various animal species (e.g., rodents), which are usually used in biological experiments. Thus a parametric identification is a challenging problem for this class of models.

Results
In this section we present some time series plots obtained by computer simulation.
For simulations the MATLAB program for modeling delay differential equations DDE23 was applied 52 . Its event location tool was used to model episodical reset of V(t). DDE23 uses the the methods of steps. Identical results were obtained when the integration was performed using a 4th-order Runge Kutta algorithm (XPPAUT 53 ).

Simulations for an exponential function F v (V).
Consider the hypothalamic input defined by (9). Let the shape function F v (V ) = e −k s V with k s = 0.5 . The circadian input is defined by (8) with t 0 = 30 min. The interval between firings varies from 51 min (at peak) to 74 min (at nadir). For similar experimentally obtained figures see 16 .
Notice that the hypothalamic oscillator and the ultradian pituitary-adrenal oscillator are coupled in both amplitude and frequency (however the maximum in CORT is phase shifted from that in ACTH). The coupling strength depends on the coefficients k ah , k vh . The pattern of oscillations for nominal parameters is shown in Fig. 2. It is seen that ACTH and CORT profiles combine at least three rhythms with the periods 20 min, ∼ 60 min and 1440 min (cf. 35 ).
The parameters were chosen to obtain the inter-burst interval close to the value of a physiological ultradian rhythm (that ranges from one hour to two). Note that the fluctuations in CORT are more regular than those observed experimentally (see e.g. 58 ). This is a consequence of the simplicity of the integrate-and-fire model we use for pituitary pulsativity. The basic ultradian and the CRH rhythms are superimposed. This may result in an additional complexity of the oscillations spectrum, since ACTH can be secreted not in spikes, but in bursts. Simulations show that this complexity reduces when the ultradian period is a multiple of the CRH period T r . (In this case the ACTH oscillation will be close to periodic with one impulse in the least period.) In other cases ACTH will have multiple impulses in the least period, i.e. bursting is observed. This is consistent with what is known about hormonal secretion in the pituitary 13,43,44,54 . With the increase of the threshold , frequency of the ultradian oscillation decreases. If the threshold is large enough, it is never reached and the system comes to an equilibrium. Note that the maximum in CORT is shifted by 15 min from the maximum in ACTH. (V). Now assume that F v (V ) = V e −k s V with the decay rate k s = 0.1 . This case will be used to compare our simulations with some experimental profiles obtained in 45 for adult male Sprague-Dawley rats. Simulated plots are shown with thick black lines, while experimental results are demonstrated with thin red lines.

Simulations for an alpha function F v
In Fig. 3A we plot a CORT profile for the circadian hypothalamic drive H(t) = H r (t)H c (t) with the circadian peak shifted to t 0 = −80 min. The simulated interval between firings varies from 71 min (at peak) to 87 min (at nadir). The comparison is made with the mean CORT concentration for the cohort of control (intact) rats (in read).
For Fig. 3B we simplify the dynamics and assume that we can neglect the modulation from the circadian SCN clock. Then the hypothalamic input can be defined not by (9), but as H(t) = H r (t) , where H r (t) is given by (7). Our modeling result is compared with Figure 3C 45 experimentally obtained for rats with SCN lesions. Our plot shows more regularity, but generally follows the experimental pattern. The basic ultradian signal has the interpulse period of approximately 76 min.

Ultradian oscillations in an isolated pituitary. Consider the case when the pituitary gland is isolated
from the hypothalamus and the adrenal gland. Here we will model the shape in vitro experiment described in 27 , where seven fetal and two adult human pituitaries were placed in a flow-through perfusion chamber and oxygenated. It was found in 27 that ACTH still exhibits a pulsatile ACTH release with low amplitudes and high frequencies. Assume that the inequality is satisfied, i.e. the firing threshold is sufficiently small, and H ≡ 0 , C ≡ 0 . Here it will shown analytically that system (1)-(5) has a globally stable periodic solution. If (10) does not hold, the system solutions approach an (10) �k v < V 0 Thus t n = t 0 + nT , n ≥ 0 , and the function V(t) is T-periodic for t ≥ t 0 .
In the case k v → +0 inequality (10) is satisfied, so V(t) is periodic. It follows a saw-tooth pattern where T = �/V 0 . Assume that (10) is valid and consider Eq. (11). Then The discrete-time mapping A(t n ) → A(t n+1 ) defined by (16) has a fixed point and for any initial value A(t 0 ) we have A(t n ) → A * as n → +∞ . It is easily seen that if we take the initial value A(t 0 ) = A * , then A(t) defined by (15) is T-periodic. www.nature.com/scientificreports/ If (10) does not hold, i.e. k v ≥ V 0 , then V (t) < � for all t ≥ t 0 and the threshold is never reached. Hence V(t) has no jumps, it is monotonously increasing for t ≥ t 0 and i.e. the solution tends to an equilibrium.

Discussion
Our mathematical model captures the main features of the HPA hormonal system, including ultradian and circadian rhythms. Our investigations were mainly inspired by pioneering works 16,42 and subsequent publications developing this approach. However, unlike these studies, we assume that the CRH drive is not constant, but periodic, with the period less than the ultradian one. Notice that unlike 16 , in our model the time delay is not critical for oscillatory behavior but only introduces a small phase shift in the secretory activity of the corticotrophs. In particular, the delay-free system is also oscillatory. Oscillatory behavior can be lost only in the case of too large values of the threshold (see Supplementary Material for exact estimates).
In our model the pulsatile exocytosis of ACTH is related to the intrinsic excitability of corticotrophs. Although corticotrophs can generate isolated spikes and bursts, it is the prolonged depolarizations associated with pseudoplateau bursting that is associated with ACTH exocytosis 13,55 . Large conductance calcium-and voltage-activated potassium (BK) channels are important for this type of bursting 13,56 . An important observation is that, at least in vitro, > 90 % of corticotrophs exhibit only spontaneous spiking activity 55 . Thus it is possible that exocytosis in vivo is related to coordinated dynamic activity in functional networks of corticotrophs 57 . Here we have assumed that bursting in the ACTH release can be explained by imposing a high-frequency CRH rhythm on the intrinsic ultradian rhythm of the ACTH. In this way our model realistically simulates the interplay between the rhythms of the PVN and the anterior pituitary.
In the model we hypothesized the existence of the three oscillation clues related to the hypothalamus and the anterior pituitary. It seems that the dynamics of a hormonal systems depends on the strength of coupling between them. It is assumed that for the HPA axis the coupling between PVN and pituitary (determined by the parameters k ah , k vh ) is rather weak. Thus the intrinsic pulse generator located in the anterior pituitary plays the leading role and determines the main hormonal rhythm. The system's dynamics is also significantly affected by the drive from the pituitary pulse generator to the adrenal gland (described with the coefficient k ca ).
Here we combine mathematical tools that were traditionally used in the modeling of endocrine systems, such as ordinary differential equations and delay differential equations, with the mathematical technique of impulsive differential equations used in describing neural activity. Such a synthesis enables explaining impulsive signalling from the brain to endocrine glands in an understandable and convenient way.
The pulses generated by the HPA axis are more irregular than those generated by our simple model (compare, for example, 12,45,58 ). Thus, the mechanism for pulse generation must be more complex than that of a simple integrate-and-fire neuron. Indeed, in vivo imaging studies suggest that hormone release by the pituitary gland is coordinated with the microcirculation 59 . In any case, it is clear that the interesting mathematics in the HPA axis occur at the level of the pituitary gland.
Our model describes the function of the HPA axis on a time scale of 24 hours. However, it is known that in conditions of chronic stress and certain psychiatric disorders the HPA axis can become dysregulated over times scales of weeks to months. Recently it has been suggested that this dysregulation arises because of slow changes in the functional masses of the corticotrophs and the adrenal cells, and a relevant mathematical model has been proposed 60 . Applying this effect to our model implies that the parameters k ah , k vh , k ca can vary in time and follow their own dynamics. Our future research will be directed towards exploring the effects of chronic and acute stress in the context of the simple model we have developed in this study.

Code availability
The code of our MATLAB program is available on request.